DIFFUSION OF A FLUID THROUGH A VISCOELASTIC SOLID 
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Abstract. This paper is concerned with the diffusion of a fluid through a viscoelastic 
solid undergoing large deformations. Using ideas from the classical theory of mixtures and a 
thermodynamic framework based on the notion of maximization of the rate of entropy pro- 
duction, the constitutive relations for a mixture of a viscoelastic solid and a fluid (specifically 
Newtonian fluid) are derived. By prescribing forms for the specific Helmholtz potential and 
the rate of dissipation, we derive the relations for the partial stress in the solid, the partial 
stress in the fluid, the interaction force between the solid and the fluid, and the evolution 
equation of the natural configuration of the solid. We also use the assumption that the 
volume of the mixture is equal to the sum of the volumes of the two constituents in their 
natural state as a constraint. Results from the developed model are shown to be in good 
agreement with the experimental data for the diffusion of various solvents through high 
temperature polyimides that are used in the aircraft industry. The swelling of a viscoelastic 
solid under the application of an external force is also studied. 



1. Introduction 

Several materials in the areas of of polymer mechanics, asphalt mechanics, and biome- 
chanics that show non-linear viscoelastic behavior swell in the presence of a fluid. Fo r 



i nstan c e, polyimides which are known to s how viscoelastic solid-like response (see iBhargava 
( 120071 ) . iFalcone and Ruggles-Wrennl ( 120091 )) are used in the aerospace ind ustry due to their 
good performance at high temperatures (also see iGhosh and Mittal these materials 

in their service environment are known to swell in the presence of moisture. In addition, 
asphalt based materials (that are well kno w to show non-lin ear viscoelastic fluid-like behav- 
ior) degrade in the presence of moisture ( Kim et aO ( 2004 )). Diffusion of biological fluids 
through biological materials is another application wherein typically nutrition is provided 
by the fluid that diffuses, and t he amount of the stress or strain in the solid can control 
the chemicals that are released ( Rajagop"al ( 20071 )). Thus, there is a considerable interest 
to understand how such viscoelastic materials deform and swell d ue to d i ffusion of a fluid. 
Study of such a phenomeno n is also of interest in geomechanics (Cohen ( 19921 )) and food 
industry dSingh et all t004 )). 



It is well known that the Darcy's and Fick's equations (iDarcyi ( 118561 ); iFick! ( 118551 )) that 
are extensively used cannot predict swelling of the solid as well as the stresses in the solid. 
In fact, Darcy's equation is an approximation of the balance of linear momentum of the fluid 
going through a rigid solid. To capture the swell ing phenomena, s everal works have been 
done using mixture theory (see review article by iRajagopall ( 120031 )) and using variational 
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principles ( jBaek and Srini vasal ( 12004] ) ). These models have been shown to match well with 
experimental data for swelling of rubber (that show elastic response) due to the diffusion of 
various organic solvents. 

In the area of diff usion of a fluid throug h viscoelastic materials, some of the earliest works 
were by Biotl ( 1956 ) and Weitsman ( 1987 ). who have used linear viscoelasticity. In deriving 
their models, they have used the fact that the fluxes and affinities ar e related through fin 
ear phenomenonological relations. Later on, Cohen and co-wor kers (ICohen and White Jr. 



(199 ll):lEdwards and Cohenl (Il995h l. and Durning and co-workers (jHuang and Durningl (ll997l ) 



Cairncross and Durning ( 19961 )) have also recognized the importance of studying the diffu- 



sion of solvents through polymeric materials that show viscoelastic response. They have 
coupled diffusion and viscoelasticity by adding terms to the flux of the diffusing fluid that 
depend on the stress in the solid. Such an approach does not have a thermodynamic ba- 
sis. Furthermo re, these models developed are only one-dimensional in nature. Recently, 



Liu et al.l (120051 ) have used the model developed by Cohen and co-workers to study the effect 



of various viscoelastic parameters on diffusion. They have shown that comparable relaxation 
times of polymer viscoelasticity and diffusion of a fluid results in non-Fickian behavior. 



1.1. Main contributions of this work. Our main goal in this paper is to develop a 
thermodynamic framework to model diffusion of a fluid through a viscoelastic solid and 
we shall mainly focus on t he swelling of po l yimides. We use i deas from mixture theor y 
see iTruesdell et al.l d200l. ITruesdelll (ll957ah. iTruesdelll (1l957bh. lAtkin and Crainej (|l976f) . 



Bedford and Drumhellerl (Il983l) . iGreen and Naghdil (ll969l ). ISamohvll ( 1l987f ). iBowenl (1l982h . 
Rajagopal and T ao (19 95J) for details) and irreversible thermodynamics to build such a frame- 
work. In lKarra and Rajagopall ( 120101 ). a framework that can be used to predict the non-linear 
viscoelastic response of polyimides under various temperat ure and loading condit i ons, h as 
been developed. In the current paper, such a framework in iKarra and Rajagopall (120101 ) is 
extended to incorporate diffusion of a fluid and to model the swelling phenomenon. 

The thermodynamic framework in the current work uses the notion of evolving nat- 
ural configuration that has been used to model a variety of phenomena including clas- 
sical plasticity, viscoelasticity, multi-network theory, superplasticity, twinning, etc. (see 
Rajagopal and Srinivasal (j2004al ) for details). The evolution of such a natural configuration is 
determined by maximizing the rate of entropy production (with any additional constraints). 
We const it utively prescribe forms for the Helmholtz potential of the mixture and the rate of 
dissipation (which is the product of density, temperature and the rate of entropy production) 
due to mechanical working, diffusion and heat conduction. The final constitutive relations 
are then derived by maximizing the rate of dissipation under appropriate constraints. In such 
an approach, one need not assume linear phenomenonological relations between the flux and 
the affinites, and thus our framework is more general. It has also been shown recently that 
if one chooses a quadratic form for the rate of entropy production in terms of affinities, and 
maximizes the rate of e ntropy production with re s pect to the affinities, one can arrive at the 



maximizes tne rate oi e ntropy production witn re s pect to tne ammties, one ca 
Onsager's relations (see Rajagopal and Srinivasa (j2004bl ) for further details). 



An initial boundary value problem is considered wherein a viscoelastic solid is held between 
two rigid walls, and immersed in a fluid. Using the model developed in the current paper, 
free swelling of such a solid and its swelling under the application of external force i.e., stress- 
assisted swelling are studied. The numerical results for free swelling of the viscoelastic solid 



are compared with experimental data for diffusion of different solvents through PMDA-ODA 
and HFPE-II-52 polyimide resins. 



1.2. Organization of the work. In section (F2J), the kinematics required in this paper are 
documented. In section ([3]), the constitutive assumptions are specified and the constitutive 
relations are derived. We shall also show that our constitutive relations reduces to the 
equations for diffusion through an elastic solid derived using theory of mixtures when certain 
parameters take special values. An initial boundary value problem is set up using our 
model in section The boundary conditions used, the non-dimensionalization scheme, 
and comparison of the numerical results with experimental data are given in ( 14. ip . ( 14. 2p . 
and ( 14. 3p . respectively. Final concluding remarks are given in section (jSJ). 

2. Preliminaries 

Let us consider a mixture of a Navier-Stokes fluid and a viscoelastic solid. We shall assume 
co-occupancy of the constituents, which is the central idea in theory of mixtures. It is based 
on the notion that at each point x in the mixture, at some time t, the two constituents exist 
together in a homogenized fashion, and are capable of moving relative to each other. We 
shall denote the quantities associated with the fluid through the superscript / and use the 
superscript s for the solid. Now, we shall define the motion x l f°r the i-th constituent of the 
mixture through 

x = X l (X\t), i = f,s, (2.1) 

where X 1 is the material point of the i-th constituent in its reference configuration. We 
shall assume that the mapping \ l is sufficiently smooth and invertible at each time t. The 
velocity associated with the i-th constituent is defined as 

= (2 ' 2) 

and the deformation gradient through 



- dx ■ < 2 - 3 > 

Let K t denote the current configuration of the mixture and let Kr, k t denote the reference 
configurations of the solid and the fluid respectively. Also, let K p [t) denote the natural 
configuration of the viscoelastic solid (see Fig. ([!])). Such a configuration is attained by the 
solid body upon instantaneous removal of external loading. For a Navi er-Stokes fluid, t he 
natural configuration is same as the current configuration of the fluid ( Rajagopal ( 1995 )). 



Now, if F\ i — f, s is the gradient of the motion (usually known as deformation gradient) 
X 1 (X l ,t), and if F S K is the gradient of the motion of the viscoelastic solid from K p (t) to 
K t , then 

GS= K w )~ ljrs - ( 2 - 4 ) 

The density p and the average velocity (also known as barycentric velocity) v of the 
mixture are defined by 

p = Xy, t> = ij>v. (2.5) 



reference configuration of the solid current configuration of the mixture 




reference configuration of the fluid 

Figure 1. Illustration of the various configurations of the viscoelastic solid 
and fluid components in the mixture. 



We define the following derivatives for any scalar quantity 0* by 

d<\> % dcj) 1 d l (j) 1 d(j) 1 . , ,,. dcj) 1 d<\> % dcj) 1 . , 

^F = ^ -dT = -dt> grad(0) = ^' d5c = d5c' (2 ' 6) 

where 

<j> i = P(x,t) = i i (X\t). (2.7) 

Hence, 

~df ~ ~dt 
and we shall also define the following 



+ grad(0*) • v\ (2.i 



f t =^+grad(0)... (2.9) 

The velocity gradient for the i-th component L l and the velocity gradient for the total 
mixture L are defined by 

L* = grad(^), L = grad(u). (2.10) 
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The symmetric and anti-symmetric parts for the velocity gradients L l , L are 



D l 



(Li 



D = - 

2 



W =2 



(2-11) 



The left Cauchy-Green stretch tensor B s v , B s p ^ and their principal invariants are defined as 



B 



G 



F s g (F s gY, 



B 



jps I pis 

K p{t) V K p(*) 



T 



(2.12) 



I/;; tr(B<), n BS = i{[tr(B5)] 2 -tr[(B5) 2 ]}, = det (Bj) , j = G,p(t), 

(2.13) 

where det(.) is the determinant of a second order tensor. We shall next note the balance 
laws. 

The balance of mass for the i-th constituent without any mass production is given by 

1 +div(pV) =0, 



at 



(2.14) 



where p % is the mass density of the z-th constituent and div(.) := tr(grad(.)) is the divergence 
operator with tr(.) meaning the trace of a second order tensor. The summation of Eq. (12.141) 
over i along with Eq. (12. 5p leads to 



— + div (pv) 



0. 



The balance of linear momentum for i-th constituent is 



(T*) T +p i b i + m i 



(2.15) 



(2.16) 



where m 1 is the interaction force on the i-th constituent due to the other constituents, b l 
is the external body force on the i-th. constituent, T l is the partial Cauchy stress tensor 
associated with the z-th constituent related to the surface traction on the i-th constituent t % 
through 

t i =(T i ) T n, (2.17) 
where n is the surface outward normal. From Newton's third law, we have 

^m* = 0. (2.18) 

i 

For mixtures, the balance of angular momentum, in the absence of body couples requires 
that the total Cauchy stress of the mixture be symmetric i.e., 



T = T T , where T = ^T* 



(2.19) 



although the individual partial stresses T l could be non-symmetric. Now, the balance of 
energy for the i-th constituent is given by 

t d 



P 



(It 



div (T l ■ v l - <f ) + pV + pty ■ v l + E i + m l ■ v l 
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(2.20) 



where e l , q l , r l are the specific internal energy, heat flux, radiant heating associated with 
the z-th component and E l is the energy supplied to the z-th constituent from the other 
constituents. 

Now, taking the scalar multiplication of Eq. ( I2.16P with v l and subtracting the resulting 
equation from Eq. ( I2.2(jp . we arrive at 

de l 

p l — = T V - div(g l ) + pV + E\ (2.21) 

Using e l = ip 1 + 9rf, where if} 1 , rf are the Helmholtz potential and specific entropy of the z-th 
constituent, and with 9 being the common temperature of the constituents at a point in the 
mixture, Eq. (I2.2ip along with Eq. ( I2.14p results in 

^ (pV ) + div (p^V) = \T l ■ V - div (£\ - l<f • grad(tf) + ipV + ±& 

11 +V 1 -jt)- (2-22) 



9 V dt dt 

Now, using the fact that rj 1 = —-£-, we can establish the following result: 

-dT + 11 H-^T' Wit -\-df- Wdi ) + v ' [ gTad ^ } " "^ gradW 



dt 



(2.23) 

9 fixed 



where the subscript ll 9 fixed" means that the derivative is to be taken keeping 9 fixed. We 
shall define 

<7=Xy> r = ij>V. (2.24) 

i P i 

Using the relation Eq. (I2.23P in Eq. (I2.22p and summing over i, along with Eq. (I2.24p . we 

get 

| (E PVJ + div PV«*J = \ E ^ - (I) - • grad(^) + p Q 

^E^E^(f) • (-) 

Eq. (I2.25P is the balance of entropy with the rate of entropy production £ being 

We shall assume that the total entropy production can be additively split into entropy pro- 
duction due to thermal effects i.e., conduction (£ c ), and entropy production due to internal 
working and mixing (( m ). We shall also require that each of these quantities be non- negative, 
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so that the rate of entropy production ( is non-negative and the second law of thermody- 
namics is satisfied automatically. This implies that 

q ■ grad(0) 

Cc : = - g2 > 0, (2.27a) 

We shall choose q = —k(p,9)grad(9), k > 0, so that Eq. (I2.27bl) automatically satisfies. 
Also, if we define the rate of dissipation £ m := 9( m , then 

£ Ti -^+E^-E^(^) =e m >o. (2.28) 



fixed 



Assuming 

E^^ + E^W-t^O, (2.29) 

i i 

Eq. (I2.28P can be re-written as 



) fixed 



Now, 



= p^ + div^pV> i -^), (2.31) 



where ^ := - E^ pV l is the average Helmholtz potential of the mixture. 



Finally, from Eqs. (I2.3ip and ( I2.30p . we arrive at 



+ div 



9 fixed 



(2.32) 



Assuming that all the components have the same Helmholtz potential Eq. ( 12.321) reduces 

to 

E^'^-E^-^-f^f) =U (2-33) 



fixed 



where we have used Eq. (12.51) . The second law of thermodynamics is invoked by ensuring 
£m > 0. The preliminaries discussed so far are sufficient for the derivation of the constitutive 
equations in section (J3J). 
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3. Constitutive assumptions 
We shall assume that the specific Helmholtz potential for the mixture is of the form 



ib = ib ( 9, Ins ll B s IIIds , l B s ,ll B s ,IIIrs ), 



(3.1) 



and so 



dib d s ib 

- = — + («-«)• grad(^) 



(£0 
#0 



) fixed 



dh 



p(t) 



I- 



#0 



+ 



p(*) . 



ails 



B s m + IIIj 



s \-l 



p(t) 



p(t) 



(3.2) 



d! B « G dllm 



dll B . 



B 



G 



+ (v - v s ) • (grad(^)) efixcd , 



(3.3) 



where ( ) is 
Now, 



dt 



for the sake of convenience. 



F = F G s + F s G 



F^FT 1 = K^G'iGT'iF^r + F° Kp(t) G s (Gr 1 (F° Kp{t ; / 



TS rs I jps rs I rps 

L ~ L P (t) + * Kv(t) L GK* K 



pit)' 



^d s = D; (t) + - 



F s Kp(t) L s G (F s Km ) 1 + (F s Kp(t) ) T (L s G ) T (F s Kp(t) ) 



(3.4) 



In addition, 



B p{t) ~ F Kp(t) (F s ) T + F s (F Kp(t) ) T 
L p(t) B p(t) + B p(t)( L p(t)) T \ 



(3.5) 



and similarly 



B a = LqBq + B S G (L G ) T . 



(3.6) 



Assuming that the response of the viscoelastic solid from the current configuration to its 
natural configuration is isotropic elastic, we choose K p ( t ) such that 



jps _ ys 

K p(t) K p(t) ' 



(3.7) 



where V s is the right stretch tensor in the polar decomposition of F 



K p(t) ' 



Using Eqs. Q33J, (E^J), (ESJ), (ETj) in Eq. fl£3P, we get 

' fixed 



2 



<9I 



+ 1 



B 



B 



MO 



ft*' (911^ 



B 



p(0 . 



00 \2 , TTT ^ ^ 



« 4 



+ 2 



#0 9-^ 
+ Is? 



DS 



^ ,nn2 , ttt d i> 



B S G Y + HIb; 



G d\\h 



■L S G 



J G) 



+ («-«').(grad(V)) flltod , (3.8) 
In what follows, we shall assume that the reference configurations (subscript o) of the 

constituents are same as their natural states (subscript R) and so S 1 := — = 2__ — 

1 ^ P 1 r detFY R 

i 

i — s,f, where we have used the fact that p % Q = p R . This need not be true in general. 



detF 1 



We shall also assume the volume additivity constraint that is based on the fact that the 
volume of the swollen viscoelas tic sol i d is e qual to the sum of the volumes of the unswollen 
viscoelastic solid and the fluid ( Mills! ( 1966 )). In our case this constraint is given by, 

b s = 1, 



+ 



and so Eq. (12.141) can be re-written as 



^ + div(0V)=O, 



which implies 



dt 



+ div ( <t> ivi 



div (4> f v f + s v s ) = (using Eq. (EHJ) ). 
Eq. (I3.12p can be re-written as 

(p s ti(L s ) + <p f ti(L f ) + v s ■ grad(0 s ) + v f ■ grad(0 / ) = 0. 
Again from Eq. (I3.9p . we have 

grad(0 s ) + grad(0 / ) = 0, 

and hence, using Eq. (13.141) in Eq. (13. 13ft . we arrive at 

(P s tr(L s ) + / tr(X / ) + v sJ ■ grad(0 s ) = 0, 

where v s j = v s — v*, is the velocity of the solid with respect to the fluid. 
Next, we shall assume that the rate of entropy production is of the form 

im = £m (Lq, F S p[t)) L f , 8, V s j) , 

and so Eq. (12.331) along with Eq. (12.181) reduces to 

T s ■ L s + T* • Lf - m s ■ v sJ - (M) = £ m (L' a , F s p(t) , L f , 6, v sJ ) . 

V 0,1 J 9 fixed 



(3.9) 
(3.10) 

(3.11) 
(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 
(3.17) 



Using Eq. in Eq. (l3"T7j) . we get 

T s ■ L s +T f ■ L f -m s ■ 



where 



T S G := 2p 



dip 



01 



+ 1 



B 



rris 

J' 1 P(t) 



dip 



T S G-L S G 



p(v-v s ) ■ (grad(V>)) efixcd 
Z m (L G ,F s p{t) ,Lf,e,v sJ ), (3.18) 



p(*) . 



#0 $0 
+ Is; 



dim, G dll B i 



B 



dtp 



a 



[B 



pit) 

s \2 , 



dip 

^w) 2 + III ^)5nf- 



B 



b: 



(3.19) 
(3.20) 



Eq. (I3.18P with the constraint Eq. (13 . 15[) can be written as 



T s ■ L s + T f ■ L f ~m s v 



J G) 



T S g-L s g ~ p(v - v s ) • (grad^) efixcd 



+ A (0 s tr(X s ) + <^tr(l/) + v sJ • grad(0 s )) = U (L S G , F s p(t) , L f , 6, v sJ ) , (3.21) 

where A is a Lagrange multiplier. 

We shall further assume that the rate of dissipation can be additively split into the rate of 
dissipation due to mechanical working of the viscoelastic solid, the rate of dissipation due to 
the Navier-Stokes fluid and the rate of dissipation due to diffusion of the fluid, with specific 
forms as follows: 



U (L' G , F s p{t) , Lf, 6, vj) = £ (L G , B s p{t) ,6) + vD* ■ Df + a(6)v sJ ■ v sJ . 



(3.22) 



Then, from Eq. (I3.22p and Eq. (13.211) . we arrive at 



T s ■ L s + T f ■ L f - m s ■ v 



pit) 



[L s -L G )-T G -L G -p{v-v> 



5rad('0)) e 



fixed 



+ A {<p s tr(L s ) + / tr(X / ) + v sJ • grad(0 s )) = £ {L' G , B s p{t) , 6) + vD* ■ D f + a(6)v sJ ■ v sJ , 

(3.23) 

which can be re-written as 

L s ■ [T s + \<p s I - T s p{f) ] + ^ ■ [Tf + X(p f I - uDf] + (T s p{t) - T G ) ■ L S G 

+ v sJ ■ [~m s + A grad(0 s ) - a(6)v sJ + pf (grad(^)) efixcd ] = £ (L S G , B s p(t) ,9) , (3.24) 

using the fact that p(v — v s ) = —pfv s j. Since, the right hand side of Eq. (13.241) does not 
depend on L s , Lf and v s j, we have 

T s = -\<j> s I + T s p(t) , (3.25) 

T f = -\(f)fl + isDf, (3.26) 

m s = A grad(0 s ) - a{9)v sJ + pf (grad(^)) efixed , (3.27) 

and so Eq. (13.241) reduces to 

{Tfa-Tb) -L G = i(L G ,B s p{t) ,9). (3.28) 

Next, we shall maximize the rate of dissipation £ with Eq. (13.281) . We shall maximize the 
auxiliary function 

$:=£ + /?[£-(T; (t) -Ty .U G ] . (3.29) 
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Now, 



^ft-W. -"■-)-«. (3,0) 



y G ^ u G 

Taking the scalar product of Eq. ( 13 .30 j) with L S G and using Eq. (I3.28p . we arrive at 

(g + 1) _ g 

and hence the evolution equation for the natural configuration of the solid is given by 



(3.31) 



(Trft) - T g) = di dL s - (3-32) 
• Lit 



dL s G 



'G 



3.1. Specific constitutive assumptions. We shall assume the following specific form for 
the specific Helmholtz potential of the mixture 

j, = A- + (B> + 4) (9 - e.) - I (9 - 0,f - 40\n (£) + ^ ^ (l Bh - 3) 

+ SeizM.^ - 3) + -*L- [(1 - # )ln(l - « - , (3.33) 

where [Aqo, Hgi, A%>o> A*pi are material parameters, is a reference temperature for the vis- 
coelastic solid, R is the gas constant, 6 the absolute temperature of mixture, Vq the molar 
volume of the fluid, and \ a mixing parameter for the particular solid-fl uid combination. The 
l ast te rm in Eq. (I3.33P is the term added to the specific Helmholtz in ( IKarra and Rajagopal 
(120101 ) ) to capture the swelling phenomenon in the solid. 
Now, 

dip 

= -(B* + c s 2 ) +cl(9- 9 S ) + c*ln (£j +4 + ^(l B s a - 3) + _ 3) 

H [(l-^)ln(l-0 s )- x (0 s ) 2 ]. (3.34) 



The internal energy e is given by 
e = ip + 9t] 

= a°- bw s + cue- e.) + § (e 2 - el) + ^(i Bb - 3) + jg(i fl . (t) - 3), (3.35) 

and the specific heat capacity C v is 

C * = % = c i 9 + c 2- (3-36) 
11 



From Eq. (ET35]) and Eq. <KW) 



Pit) 



+ 



Pr 

pR6 Jp Jq 



PrVo 

and from Eq. 1ET55)) and Eq. ( EZDj) 



2 Pp B p(t) + Pp [ l B; (t) - 3J 1 + // G (I B£ - 3) J 
[ln(l - s ) + s + X (0 S ) 2 ] /, 



2/z G .B G + /i p f I B . - 3 J I + ii G (l B s c -3)1 



+ P _RBJjJh [ln(1 _ 0S) + 0S + x(0S)2] J; 

where J1 G = BaLjbni, fip = = det (G s ), J* = det(F^ t) ). 

We shall further assume that the rate of dissipation £ is of the form 



t = 7 (e)D s G ■ d g , 



then Eq. (I3.25P becomes 



(3.37) 



(3.38) 



(3.39) 



-Xcj) s I + 



2ii p B p{t) + fi p ( Is- - 3 ) I + /i G (I B s g -3)1 



+ ^P[ln(l-0-) + ^ + X (^) 2 ] I, 



PrVo 



and Eq. (13 .32 j) reduces to 



2 ^[fl p B s p{t) -fi G B G ]= 1 (e)D G . 



The final constitutive equations are 

P 



T s = -X(t) s I + 



2jl p B s p(t) + fi p ( l Bkt) - 3 ) I + JL G (I B s g -3)1 



+ pRO ff G [ln(l - s ) + s + X (<^) 2 ] I, 
PrVo 

T f = -Xct) f I + uD f , 

m s = X grad0 s - a(6)v sJ + p f (grad^) efixcd , 



and 



2p 

P s 



[jl p B p{t) -fi G B G ]=y(e)D G} 



(3.40) 



(3.41) 



(3.42a) 

(3.42b) 
(3.42c) 



(3.43) 



being the evolution equation of the natural configuration of the solid. 

Notice that when p, G = and 7 — > 00, for the LHS of Eq. (13.431) to be finite, we must have 
D G = 0. This implies that G s = I, and hence B a u\ = B s and the solid is now an elastic 
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solid. In such a case, the constitutive equations Eq. (I3.42p . with additional assumption of 
v = 0, reduce to 



-\<p s I + 2p 



dip 



IB* 



dip 



dls* " dllB* I OIIb 

T f = -\p f I, 

m s = A grad0 s - a(6)v sJ + p f (grad^) efixed . 



^ 'B^ + IUbsJ^I 



dm 



(3.44a) 

(3.44b) 
(3.44c) 



These equations are same as the equations derived using theor y of mixtures for t he diffusion 
of a fluid through an elastic solid (see equations 3.15-3.17 in (IHron et al.l (120021 ))). 



4. Initial boundary value problem 

Let us consider the problem of compression of the viscoelastic solid body inside rigid walls 
as shown in Fig. (T5]). Let us assume that the motion of the swollen solid body is given by 



x = X, y = Y, z = f(Z,t). 
In this case, the deformation gradient of the solid {F s ) is given by 

F s = diag{l,l,p}, 

df 

where p := Let us assume a form for G s as follows 

G s = diag{l,l,s(Z,t)}, 
and the velocity of the fluid be of the form 

v' = (0,0,v(Z,t)). 

Then, 

2' 



S ; (t) =diag<M, 



(4.1) 
(4.2) 

(4.3) 
(4.4) 

(4.5) 



' l; = 2 + I - , . 



i Bk = 2 + g 2 , j; = ^, J s G = g. 



The balance of mass for the solid gives 



r det(F s ) p 
and hence from volume additivity constraint 



P f = Pi ( 1 



P 



(4.6) 



(4.7) 



(4.f 



13 



Also, we note the following relations 



Pr P \Pr 

The balance of mass for the fluid reduces to 



4 = 1 + 4(P-1), (4.9a) 

P Pr 

'= 1 + ife-lY (4.9b) 



fluid --—J ^ F(t) 

Porous filter 



.L.L.U. 



Z = H 

Swollen viscoelastic solid 



Z = 



Z = -H 

■-f-f-f-f- 

Fit) 

FIGURE 2. Schematic of the initial boundary value problem 



dp v dp , . dv 

i + ?9f + (p - 1) az=°- 

Setting v — 0, ^-component of the total stress tensor for the mixture reduces to 



J 1 — T^S _|_ rpf 

± zz — zz * zz 



A. 



= -A+(l + ^(p- I.) 



2/ip Q + /i P ^ - 1^ + hg (g 2 - 1) 



R6p 

^0 



p 



p p* 



(4.10) 



(4.11) 



(4.12) 



The balance of linear momentum for the solid and the fluid after assuming zero body forces 
reduce to 



.,d 2 f dT 2 



dt 2 dz 
dv dTl 



"'at 



dz 

14 



+ m s z , 



mi. 



(4.13a) 
(4.13b) 



Now, adding Eqs. (I4.13a[) . ( 14.13bj) . we arrive at the balance of linear momentum for the 
mixture 



P 



dt s 
dt s 



pf— = —(T s +T f ) 
1 dt dz K zz zzl 

dv dX dT s J 



dv d 



dt 



dz dz 



(4.14) 



where 



1 + 



2fl P (j) + P P - lj + fto {g 2 - 1) 



1 + 



i (a 

P \PR 



- 1 



R6p 



In ( 1-- J +~ + x\ 

pip p z 



Now, Eq. (14.13bp along with volume additivity constraint reduces to 



■ dv 



■dX 



P O = ~4 > ~r\ 1" a ~r\ V ) ~ P \ ~r\ 



dt 



dz 



df 



Of 



dip 



dz 



(4.15) 



(4.16) 



where we have used the fact that the velocity of the solid is (0, 0, 



fixed 
dt' 



Multiplying Eq. (14. 14ft with <p* and subtracting Eq. (14.161) from the resulting equation we 



get 



<P f p s 



d 2 f 



dt 2 

which reduces to 

4> f p 



+ 



■ dv 



— 1)// — = —a I — — v 1 + 



dt 2 



+ 



dt 



w dt 



df 



dt 



mil f [(hp 



dz 



> fixed 



(4.17) 



(4.18) 



where 



Pr 



(9 2 " 1) + 



Ptip P 



Pr \9 



R6p 
Wo 



1 - 



P 



In 1 



P 



1 



(4.19) 



Now, assuming that the velocity and acceleration of the solid are small compared to that of 

df d 2 f 

the fluid, we shall drop — and — — in Eq. (14. 181) . we get 

dt 



dt 



U f - l)p f ^- = av + (j) f ^r^ + P f 
dt 



■dTg 

dz 



dip 
dz 



(4.20) 



Next, we shall also assume that the acceleration of the fluid is also small and we shall drop 

dv 

— term in Eq. ( 14.201) . to get 

at 



1 

a 



dT s f 



dz 



dip 



(4.21) 
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Using Eq. (14.2 ip in Eq. (I4.10p . we arrive at 



dp 1 dp 
dt ap 2 dZ 

Also, note that 



Lf dT s J z f I dijj 



+ 



p — 1 d 
a dZ 



drit of dij 



p dZ p \ dZ 



D G = L G = diag ^ 0, 



ldg 
gdt 



and so the evolution equation of the natural configuration reduces to 

2 



7-| = 2(l + 4(p-D 

gdt \ p s r 



Peg 



(4.22) 



(4.23) 



(4.24) 



4.1. Boundary conditions. Applying boundary conditions in an initial boundary value 
problem has been an issue in mixture theory due to its basic assumption of co-occupancy. For 
instance, if traction is applied on the boundary, a natural question is how is the traction to be 
split between the solid and the fluid. To this end, the method o f spitting the traction bas ed 
on the vo lume fraction of the solid and the fluid was proposed ( Rajagopal and Taol ( 1995 )). 
Later on, iBaek and Srinivasal ( 120041 ) derived the relations for swelling of an elastic body based 
on variational principles and the boundary conditions were derived naturally. However, this 
approach assumes that the swelling is slo w, and that the relative veloci ty between the solid 
and the diffusing fluid is small. Recently, Prasad and Rajagopall ( 2006 ) have compared the 
solutions of diffusion of a fluid through a elastic slab using various boundary conditions like 
saturation boundary condition, traction splitting boundary condition, the natural boundary 
condition derived by Baek and Srinivasa, and the condition that the chemical potential is 
continuous across the boundary. Interestingly, they show that the results are insensitive to 
these different forms of boundary conditions. 

For our problem, let F(t) be the compressive force applied on the solid at Z — ±H as 
shown in Fig. (T5]) and let P^ be the pressure in the fluid at the boundaries Z = ±H, then 
we shall apply the following boundary conditions: 



-F(t) 



TL = -<P f Poo, Z = ±H, 



(4.25a) 
(4.25b) 



that is, we are assuming that the external force is borne by the solid only, while the fluid 
pressure is borne by both the solid and the fluid, and this pressure is split proportional to 
the volume fraction of the constituents. Based on these assumptions, Eq. (I4.25b[) reduces to 



A = Poo, Z = ±H. 



Eq. (I4.25ap and Eq. (|4T26|) reduce to 

-F{t)=T s J z1 Z 



±H. 



(4.26) 



(4.27) 



Note that F(t) is zero under free-swelling. 
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4.2. Non-dimensionalization. We shall use the following non-dimensionalization scheme: 



t* 



(4.28) 



where T, L are characteristic time and length respectively. If we pick fi = and define the 



non-dimensionalization quantities 0\ := ^f, 02 := L Jfy? then Eqs. ( I4.22p . (I4.24p become 

PR 



A 



+ (p* " 1) 



dp* 
p* J dZ* 

9 '1 



+ 



dZ* 



dZ* 



1 \ 1 / ffTf/* 



p" j p 



dZ* 



+ 



dZ* 



7 



1^ 



2|1 + S (P*-D 



— * / *\2 



where 

rpsf* 



0: 



(p* 



(P 



*\2 



(5 



vA2 



1 +^((^) 2 "1) 



(4.29a) 
(4.29b) 



P*/ P* 



1 



(4.30a) 



ijj* 



0i 



m* - 1) + 



px / (p 



„*\2 



A V(<?*) 



p 



The boundary conditions Eq. (I4.27P reduce to 

-F*(t*) = T S J Z \ Z* = ±H* 



(p*) 2 
(4.30b) 



(4.31) 



with F* — — , — - ss -. If we pick H as the characteristic length L, then H* = 1. The 
coupled equations Eqs. (I4.29a[) . ( 14.29b|) are solved using the staggered scheme shown in 
algorithm 1. 

The ratio of the mass of the swollen solid to its original unswollen mass can be calculated 
as follows: 



rn 



JpdV 



m J p s R dV 

f 2_ \ — dz 

J z=— 1 „s 

__ Pr 

dz 



lz=-l 

»z=l 



201 



i J z =-i L 



l + \(0i-l) 
p* 



dz. 



(4.32) 



Once, the value of p*(Z*,t*) is evaluated on the domain at various times, Eq. ( 14. 32ft is 
integrated numerically to get the mass ratio. 
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Algorithm 1 A staggered procedure for solving the coupled equations 



1: Input: ft, ft, Xi F* i P£a Ppi A*G' 7*i Time of integration tj; Time step At*; No. of 

divisions along Z* direction JV; TOLERANCE. 
2: Output: p*. 
3: Set p*° = 1, g*° = 1. 
4: while t < tf do 
5: t* = t* + At*. 
6: while true do 

7: Using p^W, flf^W, A, At, ft, ft, x, 7*, Solve for p*^ +1 ) using 

Eq. f!4.29ap in MATLAB's pdepe solver with the p*^ +1 ) at the boundaries (Z* = ±1) 
obtained by solving the non-linear algebraic equation in Eq. (I4.3ip . This non-linear 
algebraic equation is solved using the f solve solver in MATLAB. 
8: Using p* l ( l+1 \ g* 1 ^ and At, ft, ft, X-, F*, P^, p* p , pQ, 7*, Solve for g* 1 ® using 

Eq. (I4.29bp in MATLAB's ode45 solver. 
9: if ||p* i ('+ 1 ) - p*^)|| 2 < TOLERANCE then 
10: Return. 
11: end if 
12: end while 
13: p* i+l 4- p", g* i+l <- g*\ 
14: end while 



4.3. Comparison with experimental data. Fig. ([3]) shows comparison of the numeri- 
cal results to the experimental data for the ratio of swollen to unswollen (see Eq. (I4.32p ) 
PMDA-ODA (poly(N, N'- bisphenoxyphenylpyromellitimide)) due to diffusion of the sol- 
vents DMSO (dimethylsulfoxide) and NMP (N-methyl-2-pyrollidinone). In case of DMSO 
diffusing through PMDA-ODA the f ollowing material parameters were chosen: Density of 
DMSO w as chosen to be 1 096 gl ee flYang et all (j2008h ) and density of PMDA-ODA to be 
1.42 g/cc dSrinivasan et all (119941 ) and so ft = 1.3. Also, x = 0.425, (3 2 = 0.018, p* = 0.1, 
Pq = 0.1, 7* = 20 were chosen. The characteristic time (T) chosen was 10500 mm. For 
the diffusion of NMP, th e mat erial parameters chosen were: density of NMP is taken to 
be 1.02 g/cc flYang et all f l2008h ) and so ft = 1.4. Next, x = 0.6, ft = 0.016, p* = 0.1, 
Jj* G = 0.1, 7* = 20 and the characteristic time chosen was 2 45 min. The numerica l resul ts 
show good agreement with the experimental data taken from (jGattiglia and Russelll ( 119891 )). 

Next, we shall consider the diffusion of water through HPFE-II-52. The material parame- 
ters were assumed to be: x = 0-425, ft = 1.3, ft = 0.018, p* = 0.1, p* G = 0.1, 7* = 20. The 
characteristic time chose n was 2800 s. Ey en in this case the numerical results and experi- 
mental data taken from (jBhargaval (120071 )) match well (see Fig. (j4])). In all the numerical 
calculations TOLERANCE was chosen to be 10~ 4 . 

We shall now consider the problem of compression of the viscoelastic solid and study its 
effects on swelling due to diffusion of a fluid. In this numerical experiment, the solid is 
allowed to swell freely first till it saturates with fluid (upto t* = 0.5). Then, the swollen solid 
is subjected to constant compressive force of F* = 1 is applied for a time period of t* = 0.5 
and then the load is removed, and the solid is allowed to swell freely again for another time 
period of t* = 0.5. Fig. (15a|) shows that the volume of the solid gradually increases with time 
and then reaches a steady state where the volume of the solid is same everywhere and there 

18 




2000 3000 

t(min) 

(a) DMSO 



5000 



1.25 



1.15 




1.05 



40 60 80 100 

t(min) 

(b) NMP 

Figure 3. Comparis o n of the model with the experimental data from 
dGattidia and Russelll (11989^ for the diffusion of DMSO and NMP through 
PMDA-ODA (imidized at 300°C) under free-swelling condition. The charac- 
teristic times chosen were 10500 min and 245 min for DMSO and NMP, respec- 
tively. Here, 301 spatial points were used for the calculations, non-dimensional 
time step chosen is At* = 0.025. 



is no further swelling. Also, the volume of the solid near the boundary increases faster than 
that of the inner solid. Upon application of a constant compressive load in Fig. (l5bl) . the 
fluid diffuses out of the swollen solid and the volume of the solid gradually decreases until 
the volume of the solid is same everywhere. Next, upon removal of the compressive load in 
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Figure 4 . Com parison of the model with the experimental data from 
flBhargaval fl2007l )) (pg. 27) for the diffusion of water through HFPE-II-52 
under free-swelling condition. The characteristic time chosen was 2800 s. The 
parameters chosen are \ — 0.425, (3i = 1.3, (3 2 = 0.018, jit = 0.1, fi* G = 0.1, 
7* = 20. Here, 301 spatial points were used for the calculations, At* = 0.025. 

Tflytj — ^T'O 

The normalized mass is defined by , where m is the mass of the 

moo - m 

dry solid, is the steady state mass of the swollen solid, m(t) is the mass 
of the swollen solid at a given time t. 



Fig. (|6]), the solid absorbs the fluid and swells freely back to its original swollen saturation 
state. 



5. Conclusions 

We developed a systematic framework with a thermodynamic basis to develop constitutive 
relations for the diffusion of a fluid through a viscoelastic solid. A model was also derived 
using this framework by choosing specific forms for the Helmholtz potential and the rate of 
dissipation, and by maximizing the rate of dissipation. An initial boundary value problem 
was solved where we considered free swelling and swelling under the application of external 
force. We also showed that the model fits well with the experimental data for free swelling 
of polyimides. Furthermore, our work in this paper can be easily extended to study the 
diffusion of fluids in biological materials as well as in studying moisture-induced damage in 
asphalt mixes and other geomaterials that show viscoelastic behavior. 

Finally, here we shall summarize the assumptions made in this paper: 

(1) the specific Helmholtz potential of the constituents is the same, 

(2) the temperature of the constituents is the same, 

(3) the specific Helmholtz potential of the mixture depends on the temperature of the 
mixture, and the deformation of the solid, 
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Figure 5. Ratio of volume of the swollen solid to volume of unswollen solid 
(J = det(F s )) as a function of time for (a) free swelling and (b) compressive 
force F* = 1 is applied after the swollen solid reaches a saturated state due 
to free swelling. The parameters chosen are x — 0.425, f3\ = 1.3, /?2 = 0.018, 
fi* = 0.1, JIq = 0.1, 7* = 20. Here, 301 spatial points were used for the 
calculations, non-dimensional time step chosen is At* = 0.025. 



(4) the volume of the mixture is sum of the volumes of the constituents in their natural 

state, 

(5) the response of the solid from the current configuration to its natural configuration 
is isotropic and elastic, 
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Figure 6. Ratio of volume of the swollen solid to volume of unswollen solid 
(J = det(.F s )) as a function of time for free- swelling after the compressive 
force is removed. The parameters chosen are x = 0.425, fi\ = 1.3, /?2 = 0.018, 
/ip = 0.1, Hq = 0.1, 7* = 20. Here, 301 spatial points were used for the 
calculations, non-dimensional time step chosen is At* = 0.025. 



(6) the reference configurations of the constituents are same as their natural states, 

(7) the rate of dissipation of the mixture is assumed to be the sum of the rates of 
dissipation due to mechanical working of the viscoelastic solid, due to the fluid (i.e, 
due to the friction between the layers of the fluid), and due to the drag between the 
solid and the fluid. 

The following additional assumptions are made to solve the problem of compression: 

(1) the viscosity of the fluid is zero i.e., we are assuming that the dissipation due to the 
friction between the layers of the fluid is much smaller than that due to the drag 
between the solid and the fluid, 

(2) the velocity and acceleration of the solid are small compared to that of the fluid, 

(3) acceleration of the fluid is small, 

(4) the external loading is applied on the solid only, whereas the fluid pressure at the 
boundary is borne by both the solid and the fluid. 
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Appendix A. Convergence of numerical results 

Since the analytical solution for the problem is unknown, we perform an engineering 
convergence study of the solution using the described algorithm (see figure ([7])). In this 
study, the error is calculated by taking the difference between solution of various grid sizes 
(5, 15, 25, ... , upto 351 points) and the solution found using a very fine grid of 401 points at 
the point Z* = and at a time of t* = 0.5. Note that the error is propotional to logarithm 
of the spatial increment and hence the convergence rate is slow. The aim of our current work 
is not to present an optimal algorithm but to solve the coupled partial differential equations. 
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Figure 7. Engineering spatial convergence of the solution for p at Z* 
and at t* = 0.5. The time-step was chosen to be At* = 0.025. 



Appendix B. Derivation of the constitutive equation for the viscoelastic 

SOLID IN THE ABSENCE OF DIFFUSION 

Here we show that in the absence of diffusion the derived constitutiv e equations reduces to 
a varia nt of the three-dimensional standard linear solid model given in (IKarra and Rajagopal 
(120 101 )). Now, in the absence of diffusion of the fluid, we will have to drop the last term in 
(ET551I to get 

= A* + (B* + c5) (61 - fc) - f (61 - 6 s f - c^ln + ^~^ G1 V g, - 3) 

We shall assume that the solid is incompressible in the absence of fluid. The constraint of 
incompressibility is given by 

tr (D s ) = tr (D a G ) = 0. (B.2) 
The reduced energy dissipation equation of the solid reduces to 

T S ■ D S ~ ( = U (B.3) 



dt 



) fixed 



In the absence of diffusion, there will be only be dissipation due to mechanical working of 
the solid, and so the rate of dissipation in this case would be 

U = 7 (#)i^ ■ D S G . (B.4) 
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Upon maximizing the rate of dissipation using ( 1B.4j) . (1B.2j) as constraints (see (IKarra and Rajagopal 
(120101 )) for further details) we arrive at 

T s =pI + 2fi p B s p{t) , (B.5a) 

T s = XI + 2fL G B s G + TjD s G . (B.5b) 

From (IB.5ap . (1B.5bj) we have 

(p - A) I + 2ji p B s p{t) - 2fi G B s G = V D S G . 
Taking the trace of flB.6p . we get 

3 (p - A) = -2/2ptr (B; (t) ) + 2/Z G tr (B G ) . 
Hence, flB.7j) in (IB.6|) gives 



2j2 p B s p(t) - 2fi G B G = - [^tr (BJ (t) ) - /i G tr (B 8 G )] + vD & 



(B.6) 
(B.7) 

(B.8) 



The final constitutive equations for the viscoelastic solid are 

T s =pI + 2fi p B s p{t) , (B.9) 

with (1B.8|) being the evolution equation fo r the natural configura t ion of the viscoelastic solid. 
This is a variant of the model derived in ( Karra and Rajagopal ( 2010l )). 
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